#include "./position.h"
#include "../sys/rand.h"
#include <math.h>


Real mds_anint( Real x )
{
	return (Real)floor( x + MDS_F(0.5) );
}

#define MDS_ENSURE_PBC 0

Real mds_pbc1( Real x, Real L )
{
#if MDS_ENSURE_PBC == 1
	const Real y  = x - L * mds_anint( x/L );
	const Real L2 = L/2;
	assert( y >= -L2 );
	assert( y <= L2  );
	return y;
#else
	return x - L * mds_anint( x/L );
#endif
	
}


void mds_pbc_xyz(  mds_vector *r, const mds_vector *box )
{
	assert( r   != NULL );
	assert( box != NULL );
	r->x = mds_pbc1( r->x, box->x );
	r->y = mds_pbc1( r->y, box->y );
	r->z = mds_pbc1( r->z, box->z );
}

void mds_pbx_xy(   mds_vector *r, const mds_vector *box )
{
	assert( r   != NULL );
	assert( box != NULL );
	r->x = mds_pbc1( r->x, box->x );
	r->y = mds_pbc1( r->y, box->y );
}

void mds_position_initialize( mds_vector_set position, size_t n, mds_pbc_proc pbc, const mds_vector *box )
{
	assert( NULL != position );
	assert( box  != NULL );
	assert( box->x > 0 );
	assert( box->y > 0 );
	assert( box->z > 0 );
	assert( pbc );
	
	assert( n > 0 );
	{
		const Real volume  = box->x * box->y * box->z;
		const Real dr_ave  = (Real)pow(volume/n, MDS_F(1.0)/3 );
		const Real dr_min  = dr_ave * MDS_F(0.5);
		const Real dr_sq   = dr_min * dr_min;
		size_t i,j;
		
		for( i=0; i < n; ++i )
		{
			mds_vector * r_i = & position[i];
			
		TRY:
			mds_vector_rand( r_i, box );
			for( j=0; j < i; ++j )
			{
				mds_vector r_ij;
				mds_vector_vec( &r_ij, r_i, &position[j] );
				if( mds_vector_normsq( &r_ij ) <= dr_sq )
					goto TRY;
			}
			
		}
		
		for( i=0; i < n; ++i )
		{
			pbc( &position[i], box );
		}
		
	}
	
}

#if defined(__INTEL_COMPILER)
#	pragma warning ( disable : 869 )
#endif

void mds_pbc_none( mds_vector *r, const mds_vector *box )
{
}



